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Abstract 


Jet noise is a major concern in the design of commercial aircraft. Studies by various 
researchers 2-5 suggest that aerodynamic noise is a major contributor to jet noise. Some 
of these studies indicate that most of the aerodynamic jet noise due to turbulent mixing 
occurs when there is a rapid variation in turbulent structure, i.e. rapidly growing or 
decaying vortices. 

The objective of this research was to simulate a compressible round jet to study the 
non-linear evolution of vortices and the resulting acoustic radiations. In particular, to 
understand the effect of turbulence structure on the noise. An ideal technique to study 
this problem is direct numerical simulations(DNS), because it provides precise control on 
the initial and boundary conditions that lead to the turbulent structures studied. It also 
provides complete 3-dimensional time dependent data. 

Since the dynamics of a temporally evolving jet are not greatly different from those 
of a spatially evolving jet, a temporal jet problem was solved, using periodicity in the 
direction of the jet axis. This enables the application of Fourier spectral methods in 
the streamwise direction. Physically this means that turbulent structures in the jet are 
repeated in successive downstream cells instead of being gradually modified downstream 
into a jet plume. 

The DNS jet simulation helps us understand the various turbulent scales and mecha- 
nisms of turbulence generation in the evolution of a compressible round jet. These accurate 
flow solutions will be used in future research to estimate near-field acoustic radiation by 
computing the total outward flux across a surface and determine how it is related to the 
evolution of the turbulent solutions. Furthermore, these simulations allow us to investigate 
the sensitivity of acoustic radiations to inlet/boundary conditions, with possible applica- 
tion to active noise suppresion. In addition, the data generated can be used to compute 
various turbulence quantities such as mean velocities, turbulent stresses, etc. which will 
aid in turbulence modeling. 

This report will be presented in two chapters. The first chapter describes some work 
on the linear stability of a supersonic round jet and the implications of this for the jet 


1 



noise problem. 

The second chapter is an extensive discussion of numerical work using the spectral 
method which we use to solve the compressible Navier-Stokes equations to study turbulent 
jet flows. The method uses Fourier expansions in the azimuthal and streamwise direction 
and a 1-D B-spline basis representation in the radial direction. The B-spline basis is locally 
supported and this ensures block diagonal matrix equations which can be solved in O(N) 
steps. This is a modification of a boundary layer code developed by Robert Moser. 

A very accurate highly resolved direct numerical simulation (DNS) of a turbulent jet 
flow is produced. 
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Chapter 1 : The Linear Stability of a Supersonic Round Jet 


Our purpose in doing this stability problem is to provide initial conditions for the 
nonlinear computations. We want to start with the most unstable disturbance since this is 
what would be found in the natural problem. This problem has been studied by Tam and 
Hu( 1989 ) however they didn’t present the information which we require. As a result of our 
computations we have found some new results which have not been previously reported. 

The basic stability equation was obtained by linearizing the equations for inviscid 
compressible flow in cylindrical coordinates, taking small perturbations from an axially 
symmetric basic state. One can reduce this to a single second order equation for the 
pressure perturbation, which is expressed in separation of variables form as 

p = p( r ) + p^e^t-ikz-ime) 


where r is the radial coordinate. The equation for p is 


d 2 p /I 2 kdujdr 

dr 2 \ r u — uk 


1 dp\ dp /(cl i — uk ) 2 

p dr / dr V a 2 




where u(r) is the prescribed axial velocity profile, p(r) and p(r) the prescribed pressure 
and density of the base state, and a 2 (= 7 RT) is the square of the sound speed. This 
equation agrees with Tam and Hu. This equation is to be solved with boundary conditions 
which require the solution to be bounded at r — 0 and to possess only outgoing waves at 
infinity. It is an eigenvalue problem for complex u> when k and m are given. The flow is 
unstable if Imu; < 0. 

The velocity profile was taken in the same form as in Tam and Hu as a “half-Gaussian” 
function which is given by 

u = Uj r < h 

u = Uj exp ^ — ln2( — - — ) 2 j r > h 
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Half - Gaussian jet profile 



h = 0.9 r 
h + b = r 



which is sketched here in order to define the parameters b and h. We define the jet radius 
to be Vj = b + h and all computations were done with b/rj = .1, hjrj = .9. We have taken 
the temperature profile to be of a form similar to the velocity 

T = Tj r < h 

T = Too + (Tj - Too) exp ( - ln2( r - - ) 2 ) r > h 
and have taken p(r) = poo , constant across the jet. Then p is related to T by 

Poo = pRT 


and 

a 2 = aloT/Too. 

Sometimes in problems like this the temperature is related to the velocity profile by the 
Crocco-Buseman formula. But this is a viscous dominated profile which results when the 
Prandtl number is unity, and it doesnt make a lot of sense for high Reynolds number flows 
such as this. In all the work done by us so far we have taken Tj = Too so the temperature 
is constant across the jet. This could be called a “warm ” jet : the stagnation temperature 
has to be high to compensate for expansion cooling. 

We have solved this eigenvalue problem by a method described in a book by Betchov 
and Criminale. We integrate inward from a large radius, using an asymptotic formula 
to identify outgoing waves, and we integrate outward from the origin. Where the two 


4 



computations meet we must have p and dp/dr continuous. This does not occur unless 
uj has the proper value. We iterate on u> using a Newton-Raphson procedure until our 
continuity requirement is attained. The biggest problem was that the procedure doesn’t 
converge unless one starts near the proper value. This means that is was necessary to 
change parameters by small increments, starting from a known result. 

Results have been obtained for jet Mach numbers of Mj = 2 and 2.5. The result for 
Mj = 2 is shown in a figure at the end of this section. Here we present curves of growth 
rate versus k for a number of values of m (the azimuthal wavenumber). What is observed is 
that the maximum growth rate for given m increases as m increases, reaching a maximum 
at m = 4 and then it decreases again as m becomes greater than 4. This means that of 
all azimuthal wavenumbers the most unstable wave is m = 4 (and also m = —4 since only 
m 2 occurs in the equation). For Mj = 2.5 the maximum is at m = ±5. 

A wave exp(— ikz — im9) is interpreted as a helical wave on the surface of the jet, 
with the wave making an angle tan -1 (m/A;) with the axis of the jet (m = 0 is an axial 
wave). Therefore there are two helical waves of definite angle, one left-handed and one 
right-handed, which are the most unstable. This is similar to results of Sandham and 
Reynolds (1990) on the stability of a compressible mixing layer. They find two oblique 
waves with maximum growth rate which make equal but opposite angles with the flow 
direction. 

There are some implications of this for our proposed nonlinear computation. Sandham 
and Reynolds (1991) in a second paper find a single oblique unstable plane wave rolls up 
into an oblique vortex. We expect that a single unstable helical wave will roll up into a 
helical vortex. Since our helical waves propagate downstream with a speed which is about 
.66 times the jet speed, we expect a similar speed for the helical vortex. A propagating 
helix will appear like a rotating helix-as a barber pole does. This is an interesting result 
because there are observed jet modes, called spinning modes, where the noise field rotates 
around the jet axis, so that to an observer on one side the sound appears to have a time 
periodicity (Wesley and Woolley, 1975). 

Sandham and Reynolds found that the angle of the most unstable oblique wave de- 
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pends on the convective Mach number in such a way that M c cosa ~ .6, so that the 
component of the Mach number perpendicular to the wave stays constant (and subsonic). 
We find approximately the same result for the helical waves, based on only two Mach 
numbers. In the Sandham-Reynolds nonlinear computations they find that because the 
normal Mach number is subsonic the oblique vortex forms without shock waves (unlike 
their strictly two-dimensional computations which do have shocks). This is possibly very 
significant for our proposed computations. It means that helical vortices could form with- 
out shocks. This is important when we use the B-spline/spectral method because any 
shocks will have to be resolved by viscosity alone which could restrict our Reynolds num- 
ber range. 


References for Chapter 1 


1 Sandham, N. D., & Reynolds, W. C., “ Compressible mixing layer: linear theory and 
direct simulation ”, AIAA Journal , 28 , 1990, pp. 618 624. 

2 Sandham, N. D., &: Reynolds, W. C., “ Three-dimensional simulations of large eddies 
in the compressible mixing layer ”, J. Fluid Mech., 224, 1991, pp. 133-158. 

3 Tam, C. K. W., and Hu, F. Q., “ On the three families of instability waves of high 
speed jets ”, J. Fluid Mech. 201 , 1989, pp. 447-483. 

4 Westley, R., and Woolley, J. H., “ Sound pressures of a choked jet oscillating in 
the spinning mode ”, Presented as Paper 75-479 at the AIAA 2nd Aero- Acoustics 
Conference, Hampton, Va. 1975. 
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Chapter 2 : Direct Numerical Simulation of Unsteady, Compressible Round Jets 


Introduction 


Jet noise is a major concern in the design of a supersonic transport (SST) aircraft. 
Studies by various researchers 2-4 lead to aerodynamic noise as a major contributor to 
jet noise. Some of these studies indicate that most of the aerodynamic jet noise due to 
turbulent mixing occurs when there is a rapid variation in turbulent structure, i.e. rapidly 
growing or decaying vortices. 

The object of this work is to obtain highly accurate flow solutions of a turbulent 
round jet. These solutions are expected to help understand the various turbulent scales 
and mechanisms of turbulence generation in the evolution of a compressible round jet. We 
hope to use these accurate flow solutions to estimate acoustic radiation in the near-field 
region. Also the data generated can be used to compute various turbulence quantities such 
as mean velocities, turbulent stresses, etc. which may aid in turbulence modeling. 

We simulate a compressible round jet by using Fourier expansions in the azimuthal 
and streamwise. direction and a 1-D B-spline basis representation in the radial direction. 
This is an efficient and accurate way to separate out the 9 and z variables, leaving partial 
differential equations (PDEs) depending on r and t only. By using a 1-D B-spline basis 
and a Galerkin approximation we can reduce this set of PDEs to ordinary differential 
equations(ODEs) in time. This is solved by a 3 rd order Runge-Kutta time marching 
scheme. The present study uses a spectral method developed by Moser et al. 7 

We consider the temporal jet problem for two reasons: one since this configuration 
allows for the application of highly accurate spectral methods and the other because the 
dynamics of the temporal jet is not greatly different from that of a spatially evolving jet. 
The spectral accuracy helps capture smaller turbulent scales. 

Governing Equations 

The compressible Navier-Stokes equations written in cylindrical coordinates in non- 
dimensional form are, 
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where, 


d _ d d _ ld_ d d 

dx 1 dr ’ dx 2 r dd ’ dx 3 dz 


<*=)>■> m k = P u k , ^ Pr = ^ 

Re, Pr are the Reynolds number and Prandtl number, <7 p is the specific heat at 
constant pressure, and $ is the viscous dissipation (see Appendix). 


B-Spline Representation and Galerkin Formulation 


The flow variables are expanded using Fourier sums in the two periodic directions, viz. 
the azimuthal (0) and the axial (z) directions. In the non-periodic or radial direction (r) 
we use 1-Dimensional B-splines as interpolating functions. 1 B-splines have local support 
and hence lead to sparse block diagonal matrices which can be efficiently stored and solved. 
Application of boundary conditions is similar in ease to a finite element method (FEM). 

B-splines of order n are piecewise polynomials of degree n having n — 1 continuous 
derivatives. Since they have a high degree of continuity derivative quantities (like vorticity) 
can be smoothly and accurately represented. 

They have 1 degree of freedom(d.o.f) per interval unlike finite element(FE) basis which 
can have as many d.o.f as the order of the polynomial. We use B-spline bases because 
higher order FE bases may resolve waves of wavelength smaller than twice the interval 
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length due to the increased d.o.f. But this is beyond the Nyquist cut-off for the interval 
size and so these increased d.o.f do not improve the resolution (i.e. the smallest accurately 
represented scale) of the solution only thing that gets better is the accuracy (convergence 
rate is increased). Higher order B-splines on the other hand not only have better accuracy 
but also have better resolution of scales per d.o.f. 

The Galerkin method using B-splines as basis functions has been used previously by 
some researchers. 5-7 By using a Galerkin approximation we can approximate the set of 
PDEs as a set of ODEs in time. 

One can use a B-spline basis to represent the desired function f(x) as, 

OO 

/(*) = a i h j( x ) ( 4 ) 

j=- OO 

where, 

bj'(x) is a n th order B-spline coefficient, aj is the value of function at knot j. 

Using a Galerkin approximation we can write, ( bk is the weight function) 

J b k f(x)dx = '^a j J bjb k dx (5) 

Also, the derivative of f(x) can be written as, 

f( x ) = b ' 3 ( x ) ( 6 ) 

j 

But here the order of the polynomial has reduced by 1 so in order to keep the degree 
of polynomial the same we approximate the derivative as, 

f{x) « g(x) => ^ajfe'(x) « h l( x ) ( 7 ) 

3 3 

Again the Galerkin approximation is written as, 

Jb k f'(x)dx^Y, a > j b'j b k dx (8) 
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Any non-linear terms are handled in a similar way, i.e. if 

h = f g 

h = 22 dj b T j ( x ) ~ '22 a i ° k * >k 

j j,k 

22 dj J bj bi dx « 22 a 3 ° k J bj b k b[ dx (9) 

j 

The matrices (terms with the integral) on the right hand side of equations (5) and 
(8) are called the mass and derivative matrices respectively. All different combinations of 
such derivatives and other terms are computed as matrices too. These are calculated only 
once using a Gaussian quadrature and stored as opposed to a regular finite element (FE) 
basis where they can be calculated on the fly when required. 

Writing the governing equations in the Galerkin form using the B-spline representa- 
tion yields a number of matrices similar to the mass and derivative matrices, non-linear 
advection terms yield matrices similar to the non-linear matrix obtained in equation (9). 

In Galerkin form the continuity equation can be written for any kg and k z as (kg and 
k z are wave numbers in the 9 and z directions respectively), 


^ 22 J R bibi rdr = S aiG j ( mr k J r 


bi bj (r b k )' b t 


j, k, l 


where, 


, d f bi bj b k bi d [ \ 

+ -K^rng k / J - rdr + — m Zk / bi b 3 b k b t rdr 

UU Jr v OZ Jr / 

° ,t) b k (r) , mg = 22 m 0k( z ^,t) b k (r ) , 

k k 

m r = 222™> Tk (z,6,t) b k (r) , m z = ^ m 2fc (z, 6>, f ) b k (r ) 


The Fourier terms (e tka e tkz ) are included in the coefficients of the variables. So 
a k (z,6,t) = <7fc(f) J2k e e ikae+ikzZ and so on. 

Since we are using 1-D B-splines all the integrals are line integrals over the radius. 
These are computed exactly using Gaussian quadratures (doing integrals exactly takes care 
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of aliasing). The derivatives in the 9 and 2 are taken by Fast Fourier Transforming (FFT) 
into wave space and multiplying by the appropriate wave numbers kg and k z , and then an 
inverse FFT is applied to bring it back to physical space. 


The momentum and energy equations can be written in a similar manner. 



Fla 2 Computat ional Domain 

Numerical Formulation 

Writing the flow equations in a manner as discussed above results in a linear system 
of coupled equations to solve simultaneously at each time step. Since these B-splines (of 
order n) have local support onn + 1 knot(node) intervals (see fig.l) we get a 2n + 1 block 
banded matrix system. 

Mf = R (11) 

where, 
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M is the resulting mass matrix, / is the column vector of nodal values we solve for, 
and R is a column matrix resulting from the RHS of the governing equations. 

Time integration is carried out by a 3 rd order Runge-Kutta scheme. 


Regularity Requirements 


In the cylindrical coordinate system the origin (r — 0) is source of concern since some 
of the functions do not remain analytic as they have a V in the denominator. From a 
mathematical point of view the flow variables should be single valued and finite. To enforce 
this the polynomial expansion functions must satisfy some regularity requirements. 8 

The z-component of the velocity should be represented as, 


u z (r ; m, k ) = a(m, k ) P z (r 2 ; m, k) e* m 0 e l kz , m = all integers, 

( 12 ) 

Pz{ 0 ; m,k ) = 1 , 

where P z (r 2 ; m, k ) is a polynomial in r 2 . 

Scalars and z-components of all vectors should be represented in a similar manner. The 
6 and r-components of the vectors are dependent on each other and should be represented 
as, 

u r (r; m,k) = b(m,k) r^ -1 P r (r 2 ; m,k ) e 1 m6 e* kz , 
uo(r ; m,k) = c(m,k) r^ -1 P$(r 2 -, m,k ) e* rn6 e l kz , 

c(m,k) = ib(m,k ) form >1, (13) 

c(m,k ) = — ib(m,k ) form < — 1, 

P r (0; m,k) = Po( 0; m,k) = 1, 

for m = 0 b(m, k) and c(m, k) are unrelated. 

Enforcing these conditions gives rise to a set of constraint equations which then replace 
some rows in the mass matrix and suitably modify the RHS vector R. 

We can obtain the constraint equations as follows, for a quadratic B-Spline and any 
scalar u , 
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u(r; m, k) = a(m, k) P(r 2 \ m , k) e* me e l kz , m — all integers 


Only non-zero derivatives allowed are \m\ + 2 j M j > 0 

For quadratic B-splines (see fig.l), at r = 0 , 3 splines have support at the origin. 
Order of non-zero derivatives for the splines is, 

Spline 1 => 0, 1, 2 
Spline 2 =>- 1, 2 
Spline 3 =>■ 2 

All other splines are zero at r = 0. 

Consider spline expansions as Yli a % h , where bi are the spline coefficients. We have 
in any interval, 


bj = 1 =» bi + b 2 + ^3 = 1 
i 

b[ + b' 2 = 0 =► b[ = -b' 2 

Now consider different values of the azimuthal wave number m, 
for m = 0 

Only non-zero derivatives allowed are 0, 2, . . . even powers, so all odd powered deriva- 
tives should be forced to zero. At r = 0 splines having non-zero first derivative are splines 
1 and 2, so, 


Olfcj -f* 02^2 — 0 


(14) 


for rri = 1 


Only non-zero derivatives allowed are 1, 3, ... odd powers, so all even powered deriva- 
tives should be forced to zero. At r = 0 splines having non-zero 0 th and 2 nd derivatives 
are splines 1 and splines 1, 2, and 3 respectively, giving constraints, 


a\ = 0, and 
02&2 + a 3^3 = 0 


(15) 
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Similarly, 
for rn = 2 


Only non-zero derivatives allowed are 2, 4, ... even powers, giving 

a\ = 0 , 02 = 0 ( 16 ) 


and 

for m = 3 

Only non-zero derivatives allowed are 3, 5, ... odd powers, so 

a i = « 2 = a 3 = 0 (17) 

or, all splines should be zero. For all other values of m we get the same result as (17). 
Now equations (14-17) are the constraint equations which will replace the appropriate 
rows in the matrix M and vector R. 

Before replacing rows which actually contain the physics of the flow care has to be 
taken to see that no information is lost. So, we need to take a linear combination of those 
rows to be replaced and add them to the remaining rows in a particular manner. For this 
the null space of the constraint equations has to be computed and the eigenvectors of this 
space are to be used for the linear combination. This is done as shown below. 

for opoZ = 7 and m = 1 

The first set of rows in the mass matrix are written as, 



n 

0 

0 

0 

0 

0 

°\ 


* 

* 

* 

* 

* 

* 

* 


0 

bl 

4 

0 

0 

0 

0 

m' = 

* 

* 

* 

* 

* 

* 

* 


0 

b\ 

4 

4 

4 

0 
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* 

* 

* 

* 

* 

* 

* 


Vo 

4 

4 

4 

4 

^6 

b f/ 


The rows with *’s are the unconstrained rows and the 6* are the bspline derivatives. 
Now to choose a set of null vectors Xj such that, 

m'Xj = 0 (18) 
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Let us choose 


0 

0 

0 \ 

1 

0 

0 

Xu 

0 

0 

0 

1 

0 

*12 

*21 

0 

0 

0 

1 

*13 

*22 

*31 / 


So we can compute the null vectors using equation 18. 

Now the mass matrix M is modified by using the null vectors for linear combinations 


row2 = row2{original) + Xu * row^ + x\2 * row 5 + £13 * row7 
row 4 = row 4 (original) 4- X21 * row 5 + £22 * row 7 
rowe = row 6 (original) + X31 * row 7 
This has to be done for all values of the azimuthal wave number m. 

This procedure can similarly be applied to all the flow vectors and scalars appearing 
in the vector /. The derivation of constraints is similar for higher order splines. 

CFL and Modal Reduction 


Another concern arising due to a cylindrical mesh is that near the origin the aspect 
ratio of the cells gets very high, so we have to reduce the number of azimuthal modes 
near the origin to maintain a good CFL number. This we will call mode suppression. So 
we effectively reduce the number of azimuthal modes to 1 near the origin and increase it 
successively such that at the outer boundary we have all the azimuthal modes. 

This reduces the accuracy but helps us increase the time-step, dt , by a significant 
amount. Another advantage of doing this is that the computational domain in Fourier 
space is significantly smaller thus allowing faster computations, i.e. For lower azimuthal 
wave number x, the computations have to carried out for more radial points, but as x 
increases fewer radial points have to be considered. Another way to look at it would be 
near the origin fewer x loops need to be considered. 
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This is to be implemented by making the upper bound of the ar-loops to be dependent 
on y, i.e. x= 1 to nx-modified[y] if the y-loop is the outer loop, or making the lower bound 
of the y-loop dependent on x if the a>loop is the outer loop, i.e. y = ny-min[:r] to ny. 


Optimizations 


The code was optimized thus, 

• Moved all major computations (rhs of governing equations) to separate subroutines 
thus reducing the load on the main program and enabling it to be unrolled and vec- 
torized more efficiently. 

• Manually unrolled the quadrature sum loops. This helped in decreasing the number 
of operations by eliminating a lot of repetitive computations due to proper regroup- 
ing. Similar and symmetric matrix multipliers were regrouped together thus reducing 
redundant calculations. 

• Vectorizing over larger loops. 

Following the above procedure has helped reduce the time taken per mode per Runge- 
Kutta step considerably. 


Boundary conditions 

Ideally, the computational domain should mimic the physical domain by including 
all free space and only having physical boundaries. But this is not possible and the 
computational domain has to be finite. So we need artificial boundaries to limit the 
computational domain. But we want these artificial boundaries to be invisible to the 
flow field so that vortices and other waves can pass through these boundaries and leave 
the domain without giving rise to spurious reflection waves. So we need non-reflecting 
boundary conditions (NRBCs). 

In solving a temporally evolving jet the inflow-outflow boundaries are made periodic 
(see fig 2). This might not be very accurate but serves the purpose of studying turbulence. 
Also it takes care of the inflow-outflow boundary conditions. This also enables us to use 
a spectral expansion in the axial direction. So the only boundaries of concern are the 
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transverse numerical boundaries. If we have good boundary conditions we can make the 
computational domain smaller and thus get a higher resolution and compute finer scales. 

Several investigators 9-14 have studied and used different types of non-reflecting 
boundary condtions. Engquist and Majda 12,13 study the wave equation and develop 
a perfectly non-reflecting boundary condition using pseudo-differential operators. But this 
equation is non-local in both space and time and is of little use for computational pur- 
poses as we would have to store data on the boundary from previous time-steps. But by 
approximating the operator they also develop a hierarchy of approximate NRBCs , with 
the increasing accuracy in passing obliquely incident waves. 

We use the first order non-reflecting boundary conditions for outflow at the transverse 
boundary, which can also be derived from physical considerations using 1-D Riemann 
Invariants. The outflow boundary condition is, 


d_ , _ d_ , 

dt P p °° c °°dt Ur 


£oo / 

" 2 R P 


(19) 


where R is radius of outer boundary, and the primed quantities represent the pertur- 
bations from mean flow and the oo quantities represent free stream or mean flow. The 
term on the right hand side is a correction for a cylindrical boundary. 

This boundary condition is supposed to ensure that there are no incoming waves from 
infinity. It works well for directly incident waves. We think that this will serve the purpose 
since we would have nearly cylindrical wavefronts incident on the boundary which is also 
cylindrical. In addition we hope to mitigate the outgoing wave amplitudes by coarsening 
the mesh gradually as we approach the boundary as was done in the 2-D preliminary work. 

Now to implement this boundary condition in our scheme we have to transform the last 
set of equations, which represent the physics at the boundary, to characteristic variables 
at each substep of the Runge-Kutta algorithm. 

This is done using the following transformation matrix from Giles. 9 
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So to implement the above b.c. we have to replace the third of the last five rows 
including the right-hand-side by the above equation since the third characteristic C3 has 
the right form for the b.c. 
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Test Case : Spherically Radiating Flow in a Cylindrical Domain 


A few simple test cases were run to validate the code and for diagnostic purposes. In 
all cases the desired steady state solutions were obtained to machine accuracy. Here we 
present a case to test the non-reflecting b.c’s. 



Fig. 3 

We consider a spherically radiating point source flow placed at the center of a cylin- 
drical domain. The solution is singular at the origin and is therefore picked up after a 
finite time and is treated as an initial value problem. This test case is very robust, since 
there exists an exact time-dependent solution with which we can compare our results. This 
test case not only tests the non-reflecting boundary conditions but also validates the entire 
code. 

The radiating bubble is considered as a perturbation to the mean flow which is, 

P = i a = 1, u T = u z = uq = 0 

The perturbation is taken as a point source solution of the linear wave equation which 
is written as a velocity potential $(t — |). Where x = \/r 2 + z 2 , c is sound speed, and t 
is the time. can be any function. 

We have used a cubic function( a gaussian function has also been tried) for <f>, 



<F = 


47TI | ) 3 (^1 - f + f ) 3 > C (^ “ ^l) < ® 

0 otherwise 


The function is plotted below. 



Fig. 4 

The fundemental perturbations are expressed in terms of the velocity potential as, 


P 


~ Poo 


dt 


o — —p 

, d$ 


Uc 


U . 


0. 

~d~z 


Plotted below are the density profiles as time progresses. The dotted profiles which 
almost overlap the solid lines are the exact solutions. From the plots we can see that 
the non-reflecting conditions work fairly well. The waves which are incident normally 
pass through the artificial boundaries with little or no reflection. In figures 5 and 7 the 
plots with t = ts show a weak reflected residual wave. This is as expected since the non- 
reflecting conditions are exact for normally incident waves, but oblique waves would cause 
some reflection. 
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A more general test case was done by shifting the origin of the point source off the 
z-axis. So the flow is not symmetric anymore and u' e is not zero. This allows more 
azimuthal modes and also the the waves reaching the boundary are not normal incident 
anymore. This also validates the flow solver. Figure 8 shows very minimal reflection for 
the off-centered test case. 



Fig. 6 
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Fig. 8 
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Preliminary Results for M = 2 Jet 


A preliminary computation with a coarse mesh (64*67*32) was done for M = 2 jet. 
The initial conditions were taken to be a “half-Gaussian” velocity profile super-imposed 
with a combination of the most unstable eigenfunctions of the stablility problem. We have 
a single helical mode with kg = 4 and k z = 3. This makes the jet length to be about 2.1 
times the jet radius. The Reynolds number based on the jet radius and velocity of sound 
for this simulation is 2500. 

Since we use instabilities in the azimuthal direction which have a wave number of 4 the 
ensuing flow will have structures which are multiples of 4. So the quarter domain problem 
can be solved, since there is quadrant symmetry. This was verified in the preliminary run 
to ensure the validity of the code. 

Below are presented the results of a quarter domain simulation on a relatively coarse 
mesh (32*67*84). So that is 32 modes in the azimuthal direction but over a quarter 
domain. The Jet unperturbed velocity is M = 2. The Re is still 2500. The simulation was 
stopped when undamped delta waves begin to grow. At present we are working on a filter 
to remove these spurious waves. 

The energy spectrum is plotted versus the axial and the azimuthal wave numbers. 
The plots show a wide spectrum of energy. The initial conditions had energy only in the 
k z = 1 and kg = 1 modes. This indicates that the energy has been cascaded to other 
modes as the jet developed. So as expected a wide range of scales can be observed. The 
solid line shown is the k ~ 5/ 3 spectrum. So we get a decade of k ~ 5 / 3 spectrum. Better 
resolution gives us better agreement with the k ~ 5 / 3 spectrum. 

Figure 11, below shows the growth in the mean thickness of the shear layer at different 
times of the simulation. We observe a four-fold increase in the shear layer thickness from 
initial to the final time at which the simulation was stopped. 

Shown in figures 12 and 13 are cross-sectional plots of the three components of ve- 
locity, vorticity and the pressure and density at the final time (t = 6.564). The formation 
of spiral structures is also visible. A flow that was initially laminar develops turbulent 
structures. A higher resolution simulation would help capture the smaller scales better. 
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Fig. 13 : r — 0 cross-sectional plots of density and vorticity components 


Shown below is a reconstruction of the jet by placing adjacent to each other flow 
solutions at different times. The phase velocity was taken to be the velocity with which 
the disturbances traverse the domain in the axial direction. So the plots adjacent are taken 
at time intervals equal to the time it takes the disturbances to travel the domain length. 
We do see some discontinuity at the zonal boundaries since they are not exactly at the 
time intervals desired. But the picture looks good enough for us to conclude that such 
a reconstruction could be used not only to see the jet structure and growth but also to 
compute the near field acoustic signature. The flow quantity plotted is the total vorticity 
magnitude. It is evident from the plot that the initially sharp shear layer grows thicker 
and has a more varied structure showing the growth and development of helical vortices. 

Plotted above the jet is the total pressure energy (p * u r ) (or the energy outflux) 
radiated through a cylindrical shell , chosen to be r = 3 in this case. This gives a rough 
estimate of when acoustic energy was radiated. From the plot we see that there is a sharp 
rise in the energy flux after t = 2. This is when the initial disturbances in the shear layer 
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propogate to r = 3. Also around t = 4 there is another fluctuation which is due the growth 
and non-linear interactions. The energy outflux increases as we proceed in time indicating 
that acoustic emissions are now due to the tubulent behavior of the jet. There is growth of 
helical vortices. We also expect to some sharp fluctuations in the plot when the breakdown 
of the helical vortices occurs. So a direct correlation can be drawn between the flowfield 
and the affected energy flux. Bearing in mind that what happens at the shear layer (r = 
1) takes another 2 units of time to reach the location r = 3 at which the energy flux is 
computed. 


Energy outflux vs time ( at r ■ 3) 



t« 0.073 t a 1.962 t = 3.834 t* 5.806 


Fig. 14 : Vorticity magnitude at different times and the total energy flux associated with this flow 
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Appendix 


The compressible Navier-Stokes equations in cylindrical coordinates 
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p is the viscosity, 
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q is the heat flux, T is the temperature, k is heat conductivity, 
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and <t> is viscous dissipation given by, 
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